setwd("/Users/jpan/Documents/GradSchool/Disseration/Data/Dibao/Replication")

################################################################################################
###### Replication code
###### Buying Inertia: Preempting Social Disorder with Selective Welfare Provision in Urban China
###### Updated: March 9, 2015
################################################################################################

################################################################################################
###### Chapter 2

overtime <- read.csv(file="chapter2.csv",head=TRUE,sep=",",check.names=F)
overtime$Date <- as.Date(overtime$Date)

o.split <- split(overtime, overtime$City)
cities <- names(o.split)
n <- length(names(o.split))

attach(o.split)

# Figure 2.1 Left panel
plot(Beijing$Date, Beijing$Dibao, xlab="Date",ylab="Dibao Line (RMB Per Month)", main ="", ylim=c(70,530), xlim=c(as.Date("1995-01-01"),as.Date("2015-01-01")),col="white")
for(i in 1:n){
	data <- o.split[[i]]
	text(as.Date("1995-01-01"), jitter(data$Dibao[1],2.2), cities[i], cex=0.7, pos=4, col=i)
	text(as.Date("2015-01-01"), jitter(data$Dibao[32],0.5), cities[i], cex=0.7, pos=2, col=i)
	points(data$Date, data$Dibao, type="l", col=i, lty=i)
}

# Figure 2.1 Right panel
plot(Beijing$Date, Beijing$Dibao_Inf2011, xlab="Date",ylab="Dibao Line (RMB Per Month Inflation Adjusted)", main ="", ylim=c(100,530), xlim=c(as.Date("1995-01-01"),as.Date("2015-01-01")),col="white")
for(i in 1:n){
	data <- o.split[[i]]
	text(as.Date("1995-01-01"), jitter(data$Dibao_Inf2011[1],2.2), cities[i], cex=0.7, pos=4, col=i)
	text(as.Date("2015-01-01"), jitter(data$Dibao_Inf2011[32],0.5), cities[i], cex=0.7, pos=2, col=i)
	points(data$Date, data$Dibao_Inf2011, type="l", col=i, lty=i)
}

# Figure 2.2
for(i in 1:length(names(o.split))){
	data <- o.split[[i]]
	plot(data$Date, data$Dibao, type="l", ylim=c(50,700),xlab="",ylab="", main=cities[i], cex.lab = 4.5, cex.axis = 4.5, xaxt="n",yaxt="n")
	axis(2, at=	c(200,600),labels=c("200","600"), cex.axis=4.5)
	points(data$Date, data$FoodExp, col="red", pch=19, cex=2)
	text(x=as.Date("1997-01-01"),y=660,data$City[1], cex=4.5, pos=4)
}

detach(o.split)


################################################################################################
###### Chapter 3

data <- read.csv("chapter3.csv", sep=",", quote="")
data$Date <- as.Date(data$Date, origin="1899/12/31") 

# Figure 3.1 Left panel
data$Year <- format(data$Date, format="%Y") 
barplot(table(data$Year),main="",ylab="Total Ad Hoc Benefits",xlab="Year", cex.axis=1.5, cex.lab=1.5,cex=1.5, ylim=c(0,500))

# Figure 3.1 Right panel
citysort <- sort(table(data$City))
par(mar=c(5, 7, 2, 2) + 0.1)
plot(as.vector(citysort),1:36,yaxt="n",xlab="Ad Hoc Benefits (2008 - 2013)",ylab="", col="black",pch=19, cex.axis=1.5, cex.lab=1.5)
axis(2,at=1:36, labels=names(citysort), las=2, cex.axis=1.1)
abline(v=mean(as.vector(citysort)), col="red", lty="dashed", lwd=2)
text(mean(as.vector(citysort)), 33, "mean", pos=4, cex=1.5, col="red")

# Figure 3.2
city.split <- split(data, data$City)
cities <- names(city.split)
attach(city.split)

for(i in 1:length(names(city.split))){
	tempdata <- city.split[[i]]
	count <- table(tempdata$Date)
	plot(as.Date(names(count)), as.vector(count), xlim=c(as.Date("2008-1-1"),as.Date("2013-6-30")),ylim=c(0,4), col="black", pch = "|", xlab="",ylab="", main=cities[i], xaxt="n",cex.lab=4.5, cex.axis=4.5, cex=2.5)
	lines(as.Date(names(count)), as.vector(count))
    text(as.Date("2008-1-1"), 0.2, cities[i], cex=4.5, pos=4)
}

detach(city.split)

# Table 3.1
rev(sort(table(data$Cat)))

# Benefits by funding source
table(data$Gov)
table(data$Gov)[3] / sum(table(data$Gov)[2:6])  # Government
(table(data$Gov)[2] + table(data$Gov)[4] + table(data$Gov)[6]) / sum(table(data$Gov)[2:6]) # SOE, GONGO, JV
table(data$Gov)[5] / sum(table(data$Gov)[2:6])  # Private

# Benefits by form
table(data$Cash)
table(data$Cash)[1] / sum(table(data$Cash))
table(data$Cash)[2] / sum(table(data$Cash))
table(data$Cash)[3] / sum(table(data$Cash))

# Benefits by limit to the number of beneficiaries
table(data$Beneficiary_scope)

# Benefits by target
table(data$Beneficiary_target)

table(data$Beneficiary_target)[6] / sum(table(data$Beneficiary_target))  # unspecified
table(data$Beneficiary_target)[1] / sum(table(data$Beneficiary_target))  # households with children
table(data$Beneficiary_target)[3] / sum(table(data$Beneficiary_target))  # elderly

# Figure 3.3 Left panel
lianghui <- read.csv("chapter3_lianghui.csv", sep=",", quote="")
  ## variable coding:
   # start: start of first meeting
   # type: whether CPPCC or NPC is first, if "either" that means the two meetings started simultaneously
   # othermtg: number of day after "start" that the second meeting begins
   # cny: date of chinese new year (DELETE?)
   # cnydiff: "start" - date of chinese new year (negative means new year prior to "start")
lianghui$start <- as.Date(lianghui$start,origin="1899/12/31")  # start of first meeting

buffer <- 10
mtglength <- 5

lianghui$pre10 <- NULL  # number of benefits 10 days before "start"
lianghui$during <- NULL # "start" to 5 days after 2nd mtg start
lianghui$post10 <- NULL  # 10 days after 2nd mtg ends
lianghui$diff10 <- NULL # difference pre - post
lianghuidiff <- matrix(data =NA,nrow=0,ncol=3) # benefit relative to lianghui (neg = pre, pos = post)

for(i in 1:length(lianghui$start)){
	city <- lianghui$city[i]
	pre_start <- lianghui$start[i] - buffer
	pre_stop <- lianghui$start[i] + lianghui$othermtg[i] + mtglength
	post_stop <- pre_stop + buffer
	pre <- data[data$City == city & data$Date >= pre_start & data$Date < lianghui$start[i],]
	during <- data[data$City==city &data$Date >= lianghui$start[i] & data$Date <= pre_stop,]
	post <- data[data$City == city & data$Date > pre_stop & data$Date <= post_stop,]

    # count of benfits before and after lianghui for each lianghui meeting
	lianghui$pre10[i] <- dim(pre)[1]
	lianghui$during[i] <- dim(during)[1]
	lianghui$post10[i] <- dim(post)[1]	
	lianghui$diff10[i] <- dim(pre)[1] + dim(during)[1] - dim(post)[1]	

	# difference in days between benefit and lianghui
    predate <- pre$Date - lianghui$start[i]
	premat <- cbind(rep(as.character(city),length(predate)), rep("pre",length(predate)),predate)
	colnames(premat) <- c("city","period","datediff")
    lianghuidiff <- rbind(lianghuidiff, premat)
    duringdate <- during$Date - pre_stop
	duringmat <- cbind(rep(as.character(city),length(duringdate)), rep("during",length(duringdate)), 
	             duringdate)
	colnames(duringmat) <- c("city","period","datediff")
    lianghuidiff <- rbind(lianghuidiff, duringmat)
	postdate <- post$Date - pre_stop
	postmat <- cbind(rep(as.character(city),length(postdate)), rep("post",length(postdate)), 
	           postdate)
	colnames(postmat) <- c("city","period","datediff")
    lianghuidiff <- rbind(lianghuidiff, postmat)
}


lianghuidiff.pre <- lianghuidiff[lianghuidiff[,2] == "pre",]
lianghuidiff.during <- lianghuidiff[lianghuidiff[,2] == "during",]
lianghuidiff.post <- lianghuidiff[lianghuidiff[,2] == "post",]

count.during <- as.numeric(lianghuidiff.during[,3])
count.during.smooth <- c(rep(-1,length(count.during)/3),rep(-2,length(count.during)/3),rep(-3,length(count.during)/3)) # note: duration of lianghui varies by localities, and average lianghui length is 3 days, smooth across 3 days for ease of plotting. Quantity of interest is not exact date of benefit distribution during lianghui
count.pre <- as.numeric(lianghuidiff.pre[,3]) - 3 # note: move pre-date 3 days to left to make room for during date on plot
count.post <- as.numeric(lianghuidiff.post[,3])

lianghuidiff.format <- c(count.pre,count.during.smooth,count.post)
par(mar=c(5, 6, 2, 2) + 0.1)
plot(density(lianghuidiff.format), xlim=c(-24,20), ylim=c(0,.1),main="",xlab="Days Before and After Lianghui Meeting", ylab="Number of Ad Hoc Benefits Distributed",cex.axis=1.5,cex.lab=1.5, lwd=2, xaxt="n")
axis(1,at=c(-24,-14,-4,0,10,20),lab=c(-20,-10,0,0,10,20), cex.axis=1.5)
abline(v=0,col="black",lwd=1.5, lty="dashed")
abline(v=-4,col="black",lwd=1.5, lty="dashed")
text(-14,0.09, "Before Lianghui", cex=1.5)
text(-1.5,0.09, "Lianghui", cex=1.5, srt=90)
text(10,0.09, "After Lianghui", cex=1.5)


# Figure 3.3 Right panel
lianghui_nocny <- lianghui[abs(lianghui$cnydiff) <= 20,]

lianghui_nocny$pre10 <- NULL
lianghui_nocny$during <- NULL
lianghui_nocny$post10 <- NULL
lianghuidiff_nocny <- matrix(data =NA,nrow=0,ncol=3)
for(i in 1:length(lianghui_nocny$start)){
	city <- lianghui_nocny$city[i]
	pre_start <- lianghui_nocny$start[i] - buffer
	pre_stop <- lianghui_nocny$start[i] + lianghui_nocny$othermtg[i] + mtglength
	post_stop <- pre_stop + buffer
	pre <- data[data$City==city & data$Date>=pre_start & data$Date < lianghui_nocny$start[i],]
	during <- data[data$City==city & data$Date>=lianghui_nocny$start[i] & data$Date<=pre_stop,]
	post <- data[data$City == city & data$Date > pre_stop & data$Date <= post_stop,]

    # count of benfits before and after lianghui for each lianghui meeting
	lianghui_nocny$pre10[i] <- dim(pre)[1]
	lianghui_nocny$during[i] <- dim(during)[1]
	lianghui_nocny$post10[i] <- dim(post)[1]	

	# difference in days between benefit and lianghui
	if(length(pre$Date) != 0){
	  predate <- pre$Date - lianghui_nocny$start[i]	
      premat <- cbind(rep(as.character(city),length(predate)), rep("pre",length(predate)),predate)
      colnames(premat) <- c("city","period","datediff")
      lianghuidiff_nocny <- rbind(lianghuidiff_nocny, premat)
	}
	if(length(during$Date) != 0){
	  duringdate <- during$Date - pre_stop
	  duringmat <- cbind(rep(as.character(city),length(duringdate)), rep("during",length(duringdate)),  duringdate)
	  colnames(duringmat) <- c("city","period","datediff")
      lianghuidiff_nocny <- rbind(lianghuidiff_nocny, duringmat)
	}
	if(length(post$Date) != 0){
	  postdate <- post$Date - pre_stop
      postmat <- cbind(rep(as.character(city),length(postdate)), rep("post",length(postdate)), 
	           postdate)
	  colnames(postmat) <- c("city","period","datediff")
      lianghuidiff_nocny <- rbind(lianghuidiff_nocny, postmat)
	}
}

lianghuidiff_nocny.pre <- lianghuidiff_nocny[lianghuidiff_nocny[,2] == "pre",]
lianghuidiff_nocny.during <- lianghuidiff_nocny[lianghuidiff_nocny[,2] == "during",]
lianghuidiff_nocny.post <- lianghuidiff_nocny[lianghuidiff_nocny[,2] == "post",]

count_nocny.during <- as.numeric(lianghuidiff_nocny.during[,3])
count_nocny.during.smooth <- c(rep(-1,length(count_nocny.during)/3),rep(-2,length(count_nocny.during)/3),rep(-3,length(count_nocny.during)/3)) # note: duration of lianghui varies by localities, and average lianghui length is 3 days, smooth across 3 days for ease of plotting. Quantity of interest is not exact date of benefit distribution during lianghui
count_nocny.pre <- as.numeric(lianghuidiff_nocny.pre[,3]) - 3 # note: move pre-date 3 days to left to make room for during date on plot
count_nocny.post <- as.numeric(lianghuidiff_nocny.post[,3])

lianghuidiff_nocny.format <- c(count_nocny.pre, count_nocny.during.smooth, count_nocny.post)
par(mar=c(5, 6, 2, 2) + 0.1)
plot(density(lianghuidiff_nocny.format), xlim=c(-24,20), ylim=c(0,.1),main="",xlab="Days Before and After Lianghui Meeting", ylab="Number of Ad Hoc Benefits Distributed\n(Excluding Lianghui close to Chinese New Year)",cex.axis=1.5,cex.lab=1.5, lwd=2, xaxt="n")
axis(1,at=c(-24,-14,-4,0,10,20),lab=c(-20,-10,0,0,10,20), cex.axis=1.5)
abline(v=0,col="black",lwd=1.5, lty="dashed")
abline(v=-4,col="black",lwd=1.5, lty="dashed")
text(-14,0.09, "Before Lianghui", cex=1.5)
text(-1.5,0.09, "Lianghui", cex=1.5, srt=90)
text(10,0.09, "After Lianghui", cex=1.5)

# Figure 3.4
res1 <- t.test(lianghui$post10,lianghui$pre10+lianghui$during)
res2 <- t.test(lianghui_nocny$post10, lianghui_nocny$pre10+lianghui_nocny$during)

par(mar=c(7, 7.5, 1, 2.5) + 0.1)
plot(1, 0.5,xlim=c(-1,1),ylim=c(0,1), ylab="",col="white", yaxt="n",xlab="Difference in Ad Hoc Benefits\n(After - Before Lianghui)", cex.axis=1.5, cex.lab=1.5, lwd=2, mgp=c(4,1,0))
axis(side=2, at = c(0.75, 0.4), labels=c("All\nmeetings", "Meetings\nstarting 20\ndays before\nor after\nChinese\nNew Year"),las=2, cex.axis=1.5)
points(res1$estimate[1] - res1$estimate[2], 0.75, pch=16, cex=2)
lines(c(res1$conf.int[1], res1$conf.int[2]), c(0.75,0.75),lwd=2)
points(res2$estimate[1] - res2$estimate[2], 0.4, pch=16, cex=2)
lines(c(res2$conf.int[1], res2$conf.int[2]), c(0.4,0.4), lwd=2)
abline(v=0, lty="dashed")
text(-0.55, 0.97, "More benefits\nbefore Lianghui", cex=1.5)
text(0.55, 0.97, "More benefits\nafter Lianghui", cex=1.5)

# Figure 3.5 Left panel
bycity <- read.csv("chapter3_bycity.csv")
  ## variable coding
   # Benefits2012: ad hoc benefits by city in 2012
   # pre10: ad hoc benefits within 10 days of start of lianghui by city (for lianghui 20 days or more away from Chinese New Year)
   # census10.pop.num.nonhan: number of ethnic groups by city, based on 2010 census
   # nonhan.conflict14: population share of minorities suspected of having propensity to protest
   # poverty: proportion of city population falling below Dibao line set by the city at end of 2004
   # fiscalrev10: fiscal revenue by city for 2010
   # census10.elder65: proportion of population over the age of 65 from 2010 census
   # Protests2012: reports of protests by city in 2012

plot(bycity$census10.pop.num.nonhan, bycity$Benefits2012, main ="",xlab="Number of Ethnic Minority Groups", ylab="Number of Ad Hoc Dibao Benefits", col="white", cex.axis=1.5, cex.lab=1.5)
text(bycity$census10.pop.num.nonhan, bycity$Benefits2012, bycity$City, cex=1.2)
lines(lowess(bycity$census10.pop.num.nonhan, bycity$Benefits2012), col="dodgerblue", lwd=2)

# Figure 3.5 Right panel
plot(bycity$census10.pop.num.nonhan, bycity$pre10, main ="",xlab="Number of Ethnic Minority Groups", ylab="Ad Hoc Dibao Benefits Prior to Lianghui", col="black", pch=19, cex.axis=1.5, cex.lab=1.5)
lines(lowess(bycity$census10.pop.num.nonhan, bycity$pre10), col="dodgerblue", lwd=2)

# Appendix Table A.1
summary(lm(bycity$Benefits2012 ~ bycity$census10.pop.num.nonhan)) 
summary(lm(bycity$Benefits2012 ~ bycity$census10.pop.num.nonhan + bycity$poverty + bycity$fiscalrev10 + bycity$census10.elder65 + bycity$Protests2012))  

# Figure 3.6 Left panel
library(Zelig)

bycity_nona <- na.omit(data.frame(cbind(bycity$Benefits2012, bycity$pre10, bycity$census10.pop.num.nonhan, bycity$nonhan.conflict14, log(bycity$nonhan.conflict14), bycity$fiscalrev10, bycity$poverty,
 bycity$census10.elder65, bycity$Protests2012)))
names(bycity_nona) <- c("Benefits2012","pre10","census10.pop.num.nonhan", "nonhan.conflict14","log.nonhan.conflict14","fiscalrev10","poverty","census10.elder65","Protests2012")

z.out <- zelig(Benefits2012 ~ census10.pop.num.nonhan + fiscalrev10 + poverty + census10.elder65 + Protests2012, model="ls", data= bycity_nona)
summary(z.out)
x1 <- seq(41,55, 0.1)
x.nonhan <- setx(z.out, census10.pop.num.nonhan = x1)
sim.nonhan <- sim(z.out, x=x.nonhan)
mu <- NULL
lower <- NULL
upper <- NULL
for(i in 1:length(summary(sim.nonhan)$stats)){
	mu[i] <- summary(sim.nonhan)$stats[[i]][[1]][1]
	lower[i] <- summary(sim.nonhan)$stats[[i]][[1]][4]
	upper[i] <- summary(sim.nonhan)$stats[[i]][[1]][5]
}
plot(x=x1,y=mu,type="l",lwd=2, xlab="Number of Ethnic Minority Groups",ylab="Expected Probability of Number of Ad Hoc Benefits", main="", col="white", ylim=c(5,21), cex.lab=1.5, cex.axis=1.5) #ylim=c(0,1)
for(i in 1:length(summary(sim.nonhan)$stats)){
	segments(x1[i],lower[i],x1[i],upper[i], lty=1)
}

# Appendix Table A.2
summary(lm(bycity$pre10 ~ bycity$census10.pop.num.nonhan))  
summary(lm(bycity$pre10 ~ bycity$census10.pop.num.nonhan + bycity$fiscalrev10 + bycity$poverty + bycity$census10.elder65 + bycity$Protests2012))

# Figure 3.6 Right panel
z.out <- zelig(pre10 ~ census10.pop.num.nonhan + fiscalrev10 + poverty + census10.elder65 + Protests2012, model="ls", data= bycity_nona)
x1 <- seq(41,55, 0.1)
x.nonhan <- setx(z.out, census10.pop.num.nonhan = x1)
sim.nonhan <- sim(z.out, x=x.nonhan)
mu <- NULL
lower <- NULL
upper <- NULL
for(i in 1:length(summary(sim.nonhan)$stats)){
	mu[i] <- summary(sim.nonhan)$stats[[i]][[1]][1]
	lower[i] <- summary(sim.nonhan)$stats[[i]][[1]][4]
	upper[i] <- summary(sim.nonhan)$stats[[i]][[1]][5]
}
plot(x=x1,y=mu,type="l",lwd=2, xlab="Number of Ethnic Minority Groups",ylab="Expected Probability of Ad Hoc Benefits before Lianghui", main="", col="white", ylim=c(-1,5), cex.axis=1.5, cex.lab=1.5) #ylim=c(0,1)
for(i in 1:length(summary(sim.nonhan)$stats)){
	segments(x1[i],lower[i],x1[i],upper[i], lty=1)
}

# Figure 3.7 Left Panel
plot(log(bycity$nonhan.conflict14), bycity$Benefits2012, col="white", xlab="", ylab="Number of Ad Hoc Dibao Benefits", xaxt="n", cex.lab=1.5, cex.axis=1.5)
mtext("Log Population Share of Minorities\nwith Higher Perceived Propensity to Protest", 1, line=2, cex=1.5)
text(log(bycity$nonhan.conflict14), bycity$Benefits2012, label=bycity$City, cex=1.1)
lines(lowess(log(bycity$nonhan.conflict14), bycity$Benefits2012), col="dodgerblue", lwd=2)

# Figure 3.7 Right Panel (VARIABLE: nonhan.conflict14, census10.pop.han.per)
par(mgp=c(1,1, 0))
plot(log(bycity$nonhan.conflict14), log(1-bycity$census10.pop.han.per), xlab="", ylab="Log Population Share of All Minorities", col="white", cex.lab=1.5,xaxt='n', yaxt='n', xlim=c(-9,0), ylim=c(-9,0)) #
mtext("Log Population Share of Minorities\nwith Higher Perceived Propensity to Protest", 1, line=2, cex=1.5)
abline(a=0,b=1, lwd=2, lty=2)
text(log(bycity$nonhan.conflict14), log(1-bycity$census10.pop.han.per), label=bycity$City, cex=1.1)
text(-4.2,-4.5, "45 degree line", cex=1.1, srt=45)

# Appendix Rable A.3
summary(lm(bycity$Benefits2012 ~ log(bycity$nonhan.conflict14)))  # pvalue 0.0616
summary(lm(bycity$Benefits2012 ~ log(bycity$nonhan.conflict14) + bycity$fiscalrev10 + bycity$poverty + bycity$census10.elder65 + bycity$Protests2012))  # pvalue 0.140

# Figure 3.8
z.out <- zelig(Benefits2012 ~ log.nonhan.conflict14 + fiscalrev10 + poverty + census10.elder65 + Protests2012, model="ls", data= bycity_nona)
x2 <- seq(min(as.numeric(names(table(log(bycity$nonhan.conflict14))))), max(as.numeric(names(table(log(bycity$nonhan.conflict14))))), 0.05)
x.nonhan <- setx(z.out, log.nonhan.conflict14 = x2)
sim.nonhan <- sim(z.out, x=x.nonhan)
mu <- NULL
lower <- NULL
upper <- NULL
for(i in 1:length(summary(sim.nonhan)$stats)){
	mu[i] <- summary(sim.nonhan)$stats[[i]][[1]][1]
	lower[i] <- summary(sim.nonhan)$stats[[i]][[1]][4]
	upper[i] <- summary(sim.nonhan)$stats[[i]][[1]][5]
}
plot(x=x2,y=mu,type="l",lwd=2, xlab="",ylab="Expected Probability of Ad Hoc Benefits", main="", col="white", cex.lab=1.5, cex.axis=1.5, ylim=c(5,20),xaxt="n") #ylim=c(0,1)
mtext("Log Population Share of Minorities\nwith Higher Perceived Propensity to Protest",1,line=2,cex=1.5)
for(i in 1:length(summary(sim.nonhan)$stats)){
	segments(x2[i],lower[i],x2[i],upper[i], lty=1)
}


################################################################################################
###### Chapter 4
###### See replication data for:
###### Sources of Authoritarian Responsiveness: A Field Experiment in China


################################################################################################
###### Chapter 5

data <- read.csv("neighborhood.csv", sep=",", quote="")

# Table 5.1  (VARIABLES: city, rc_num, resall)
# number of neighborhoods per city
table(data$city)

hangzhou <- data[data$city == "hangzhou",]
wuhan <- data[data$city == "wuhan",]
zhengzhou <- data[data$city == "zhengzhou",]
xian <- data[data$city == "xian",]

# number of RC members interviewed by city
sum(hangzhou$rc_num)
sum(wuhan$rc_num)
sum(zhengzhou$rc_num)
sum(xian$rc_num)

# number of residents interviewed by city
sum(hangzhou$resall)
sum(wuhan$resall)
sum(zhengzhou$resall)
sum(xian$resall)

# Table 5.2 (VARIABLES: bldgmgr2, name_bldgmgr2, dibao_disable, dibao_kid, dibao_jail, rent, hkrul)
table(data$bldgmgr2)      # Block captains
table(data$name_bldgmgr2) # High information extraction capacity
table(data$dibao_disable) # Dibao recipients: disabled
table(data$dibao_kid)     # Dibao recipients: HH with children
table(data$dibao_jail)    # Dibao recipients: ex prisoners
table(data$rent)  # apartments for rent
table(data$hkrul) # residents with rural residential permit

# Table 5.3 (VARIABLE: year)
table(data$year)  # Neighborhoods Age, 9999 = NA

# Appendix Table A.6 (VARIABLE: neighbors2)
library(Zelig)
summary(zelig(dibao_jail~name_bldgmgr2 + city, model="logit", data=data)) # column (1)
summary(zelig(dibao_jail~name_bldgmgr2+year+rent+hkrul+ city, model="logit", data=data)) # column (2)
summary(zelig(dibao_jail~name_bldgmgr2+neighbors2+year+rent+hkrul+city,model="logit",data=data)) # col (3)

summary(zelig(dibao_disable~name_bldgmgr2 + city, model="logit", data=data)) # column (4)
summary(zelig(dibao_disable~name_bldgmgr2+year+rent+hkrul+ city, model="logit", data=data)) # column (5)
summary(zelig(dibao_disable~name_bldgmgr2+neighbors2+year+rent+hkrul+city,model="logit",data=data)) # col(6)

# Figure 5.3 Left panel
data.jail <- na.omit(subset(data, select=c(dibao_jail, name_bldgmgr2, neighbors2,city,rent, year,hkrul)))
out.jail2 <- zelig(dibao_jail~name_bldgmgr2+year+rent+hkrul+city,model="logit",data= data.jail)
x.jail.l <- setx(out.jail2, name_bldgmgr2 =0)
x.jail.h <- setx(out.jail2, name_bldgmgr2 =1)
sim.jail <- sim(out.jail2, x= x.jail.l, x1= x.jail.h)
xl.est <- unlist(sim.jail$qi[1])
xh.est <- unlist(sim.jail$qi[2])
fd.est <- unlist(sim.jail$qi[5])

plot(c(0,0),xlim=c(0.5,3.5),ylim=c(-.375,1),col="white",xaxt="n",xlab="",ylab="Probability of Former Prisoners Among Dibao Recipients",main="",cex.lab=1.5,cex.axis=1.5)
axis(1,at=c(1,2,3),labels=c("High Info\nExtraction","Low Info\nExtraction","First\nDifference"), cex.axis=1.5,tick=F,mgp=c(3,2,0))
abline(h=0,lwd=1,lty="dashed")
points(1,mean(xh.est), pch=19, cex=1.5)
lines(c(1,1),c(quantile(xh.est, .05), quantile(xh.est, .95)), lwd=2)
points(2,mean(xl.est), pch=19, cex=1.5)
lines(c(2,2),c(quantile(xl.est, .05), quantile(xl.est, .95)), lwd=2)
points(3,mean(fd.est), pch=17, col="black", cex=1.5)
lines(c(3,3),c(quantile(fd.est, .05), quantile(fd.est, .95)), lwd=2, col="black", lty="dashed")

# Figure 5.3 Right panel
data.disable <- na.omit(subset(data, select=c(dibao_disable, name_bldgmgr2, neighbors2,city,rent, year,hkrul)))
out.disable2 <- zelig(dibao_disable~name_bldgmgr2+year+rent+hkrul+city,model="logit",data= data.disable)
x.dis.l <- setx(out.disable2, name_bldgmgr2 =0)
x.dis.h <- setx(out.disable2, name_bldgmgr2 =1)
sim.dis <- sim(out.disable2, x= x.dis.l, x1= x.dis.h)
xl.est2 <- unlist(sim.dis$qi[1])
xh.est2 <- unlist(sim.dis$qi[2])
fd.est2 <- unlist(sim.dis$qi[5])

plot(c(0,0),xlim=c(0.5,3.5),ylim=c(-.375,1),col="white",xaxt="n",xlab="",ylab="Probability of Disabled Among Dibao Recipients",main="",cex.lab=1.5,cex.axis=1.5)
axis(1,at=c(1,2,3),labels=c("State\nPenetration","No State\nPenetration","First\nDifference"), cex.axis=1.5,tick=F,mgp=c(3,2,0))
abline(h=0,lwd=1,lty="dashed")
points(1,mean(xh.est2), pch=19, cex=1.5)
lines(c(1,1),c(quantile(xh.est2, .05), quantile(xh.est2, .95)), lwd=2)
points(2,mean(xl.est2), pch=19, cex=1.5)
lines(c(2,2),c(quantile(xl.est2, .05), quantile(xl.est2, .95)), lwd=2)
points(3,mean(fd.est2), pch=17, col="black", cex=1.5)
lines(c(3,3),c(quantile(fd.est2, .05), quantile(fd.est2, .95)), lwd=2, col="black", lty="dashed")

################################################################################################
###### Chapter 6

data <- read.csv("neighborhood.csv", sep=",", quote="")

# Appendix Table A.7 (VARIABLE diffppl2)
summary(zelig(diffppl2 ~dibao_jail+dibao_disable, model="logit", data=data)) # column (1)
summary(zelig(diffppl2 ~dibao_jail+dibao_disable + city, model="logit", data=data)) # column (2)
summary(zelig(diffppl2 ~dibao_jail+dibao_disable+year+rent+hkrul+city, model="logit", data=data)) # col (3)

# Figure 6.2
data.diffppl <- na.omit(subset(data, select=c(dibao_jail, dibao_disable, diffppl2, city,rent, year,hkrul)))
out.diffppl2 <- zelig(diffppl2 ~dibao_jail+dibao_disable+city,model="logit",data= data.diffppl)
xl <- setx(out.diffppl2, dibao_jail =0)
xh <- setx(out.diffppl2, dibao_jail =1)
sim.diffppl2.jail <- sim(out.diffppl2, x= xl, x1= xh)
xl <- setx(out.diffppl2, dibao_disable =0)
xh <- setx(out.diffppl2, dibao_disable =1)
sim.diffppl2.disable <- sim(out.diffppl2, x= xl, x1= xh)

xl.jail <- unlist(sim.diffppl2.jail$qi[1])
xh.jail <- unlist(sim.diffppl2.jail$qi[2])
fd.jail <- unlist(sim.diffppl2.jail$qi[5])

xl.disable <- unlist(sim.diffppl2.disable $qi[1])
xh.disable <- unlist(sim.diffppl2.disable $qi[2])
fd.disable <- unlist(sim.diffppl2.disable $qi[5])

plot(c(0,0),xlim=c(0.5,6.5),ylim=c(-1,1),col="white",xaxt="n",xlab="",ylab="Probability of Dibao Complaints to Residential Committee",main="",cex.lab=1.5,cex.axis=1.5)
axis(1,at=c(1,2,3,4,5,6),labels=c("Former\nPrisoners","No Former\nPrisoners","First\nDifference","Disabled","No\nDisable","First\nDifference"), cex.axis=1.5,tick=F,mgp=c(3,2,0))
abline(v=3.5,lwd=1,lty="solid")
abline(h=0,lwd=1,lty="dashed")
points(1,mean(xh.jail), pch=19, cex=1.5)
lines(c(1,1),c(quantile(xh.jail, .05), quantile(xh.jail, .95)), lwd=2)
points(2,mean(xl.jail), pch=19, cex=1.5)
lines(c(2,2),c(quantile(xl.jail, .05), quantile(xl.jail, .95)), lwd=2)
points(3,mean(fd.jail), pch=17, col="black", cex=1.5)
lines(c(3,3),c(quantile(fd.jail, .05), quantile(fd.jail, .95)), lwd=2, col="black", lty="dashed")
#legend(x="topright",bty="n",legend="90% confidence intervals",cex=1.25)
points(4,mean(xh.disable), pch=19, cex=1.5)
lines(c(4,4),c(quantile(xh.disable, .05), quantile(xh.disable, .95)), lwd=2)
points(5,mean(xl.disable), pch=19, cex=1.5)
lines(c(5,5),c(quantile(xl.disable, .05), quantile(xl.disable, .95)), lwd=2)
points(6,mean(fd.disable), pch=17, col="black", cex=1.5)
lines(c(6,6),c(quantile(fd.disable, .05), quantile(fd.disable, .95)), lwd=2, col="black", lty="dashed")

# Appendix Table A.8 (VARIABLE: yingbao3)
summary(zelig(yingbao3 ~dibao_jail, model="logit", data=data)) #jail/dibao neg
summary(zelig(yingbao3 ~dibao_jail+city, model="logit", data=data)) #jail/dibao neg
summary(zelig(yingbao3 ~dibao_jail+year+rent+hkrul+city, model="logit", data=data)) #neg

